import delimited using ../input/parameter.csv, clear
mkmat _all, matrix(P)
matrix list P
clear		
	
forvalues c=1/15{
	scalar gamma = 0.003
	scalar mu_b = 0.04
	scalar mu_d = 0.01
	scalar sigma_eps = 0.5
	scalar sigma_eta = 0.01
	scalar delta_d = 1
	scalar sigma_d = 0.5
	
	scalar rho_ee = P[`c',1]
	scalar kappa = P[`c',2]	
	scalar delta_b = P[`c',3]
	scalar sigma_b = P[`c',4]
	scalar rho_bd = P[`c',5]
	
	disp rho_ee kappa delta_b sigma_b rho_bd
	
	local KK=20
	local NN=5000
	
	set obs `NN'
	gen w = (rnormal()>0)
	gen z = (rnormal()>0)
	gen e1 = rnormal()
	gen e2 = rnormal()
	gen eps = sigma_eps*e1
	gen eta = sigma_eta*( rho_ee*e1 + sqrt(1-rho_ee^2)*e2 )
	corr eps eta
	gen u1 = rnormal()
	gen u2 = rnormal()
	gen b = mu_b + delta_b*(w-0.5) + sigma_b*u1
	gen d = exp( ln(mu_d)+ delta_d*(w-0.5) + sigma_d*( rho_bd*u1 + sqrt(1-rho_bd^2)*u2 ) )
	corr b d
	sum b d, d
	
	gen vmax = - 10^8
	gen x = 0
	qui forvalues k=0/`KK'{
		gen o`k' = ( mu_b - eta )*`k' - d*z*max(`k'-12,0) - (gamma/2)*`k'^2
	//	gen o`k' = ( mu_b - eta - d*z )*`k' - (gamma/2)*`k'^2
		replace x = `k' if o`k'>vmax
		replace vmax = max(vmax,o`k')
	}
	gen y = b*x + kappa*(x>=12) + eps
	
	tab x, m
	sum x y, d	
	
	qui forvalues j=6/`KK'{
		gen x`j' = (x>=`j')
		eststo: reg x`j' x w
	}
	esttab, b(%6.3f) nose not keep(x) nostar compress
	eststo clear
	qui forvalues j=6/`KK'{
		eststo: ivregress 2sls x`j' (x=z) w
	}
	esttab, b(%6.3f) nose not keep(x) nostar compress
	eststo clear
	
	drop _all
	
	local RR=1000
	if `RR'>400{
		set matsize `RR'
	}
	matrix A = J(`RR',16,.)
	mat colnames A = x_avg x_sd ols ols_se iv iv_se gap gap_se cw cw_se tw tw_se me me_se lm lm_se
	
	qui forvalues j=1/`RR'{
		set obs 5000
		gen w = (rnormal()>0)
		gen z = (rnormal()>0)
		gen e1 = rnormal()
		gen e2 = rnormal()
		gen eps = sigma_eps*e1
		gen eta = sigma_eta*( rho_ee*e1 + sqrt(1-rho_ee^2)*e2 )
		gen u1 = rnormal()
		gen u2 = rnormal()
		gen b = mu_b + delta_b*(w-0.5) + sigma_b*u1
		gen d = exp( ln(mu_d)+ delta_d*(w-0.5) + sigma_d*( rho_bd*u1 + sqrt(1-rho_bd^2)*u2 ) )
		gen vmax = - 10^8
		gen x = 0
		forvalues k=0/`KK'{
			gen o`k' = ( mu_b - eta )*`k' - d*z*max(`k'-12,0) - (gamma/2)*`k'^2
			replace x = `k' if o`k'>vmax
			replace vmax = max(vmax,o`k')
		}
		gen y = b*x + kappa*(x>=12) + eps
	
		sum x	
		matrix A[`j',1] = r(mean)
		matrix A[`j',2] = r(sd)
		
		foreach v in x z{
			reg `v' w
			predict e_`v', resid
		}
		reg y x w, robust
		predict e_ols, resid
		scalar b_ols = _b[x]
		matrix A[`j',3] = b_ols
		matrix A[`j',4] = _se[x]
		
		ivregress 2sls y (x=z) w, robust
		predict e_iv, resid
		scalar b_iv = _b[x]
		matrix A[`j',5] = b_iv
		matrix A[`j',6] = _se[x]
			
		gen x2 = e_x^2
		sum x2
		local x2_m = r(mean)
		gen xz = e_x*e_z
		sum xz
		local xz_m = r(mean)
		gen temp = e_iv*e_z/`xz_m'-e_ols*e_x/`x2_m'
		reg temp, robust
		matrix A[`j',7] = b_iv - b_ols 	
		matrix A[`j',8] = _se[_cons]
		drop temp
		
		// First stage for OLS with IV covariate weight 
		gen xw = x*w
		reg y x xw w
		predict y_cw, xb
		predict v_cw, resid
		reg e_z x xw w
		predict z_cw, xb
		reg y_cw w
		predict yr_cw, resid
		
		// First stage for Lochner-Moretti reweighted OLS
		reg y i.x w
		predict y_lm, xb
		predict v_lm, resid
		reg e_z i.x w
		predict z_lm, xb
		
		// First stage for OLS with IV weight
		egen grp = group(x w)
		areg y, a(grp)
		predict y_tw, xbd
		predict v_tw, resid
		areg e_z, a(grp)
		predict z_tw, xbd
		
		// CW diff
		// stata refuses to run ivregress when y_cw is perfectly predicted by (x,w); manually run the regressoion 
		gen temp = yr_cw*e_z
		sum temp
		scalar b_cw = r(mean)/`xz_m'
		gen e_cw = yr_cw - b_cw*e_x
		drop temp
		matrix A[`j',9] = b_cw - b_ols
		gen temp = ( v_cw*z_cw + e_cw*e_z )/`xz_m' - e_ols*e_x/`x2_m'
		reg temp, robust
		matrix A[`j',10] = _se[_cons]
		drop temp
		
		// TW diff
		ivregress 2sls y_tw (x=z) w
		predict e_tw, resid
		scalar b_tw = _b[x]
		matrix A[`j',11] = b_tw - b_cw
		gen temp = ( ( v_tw*z_tw + e_tw*e_z )-( v_cw*z_cw + e_cw*e_z ) )/`xz_m'
		reg temp, robust
		matrix A[`j',12] = _se[_cons]
		drop temp
		
		// ME diff
		matrix A[`j',13] = b_iv - b_tw
		gen temp = ( e_iv*e_z - ( v_tw*z_tw + e_tw*e_z ) )/`xz_m'
		reg temp, robust
		matrix A[`j',14] = _se[_cons]
		drop temp	
	
		
		// LM decomposition
		ivregress 2sls y_lm (x=z) w
		predict e_lm, resid
		scalar b_lm = _b[x]
		matrix A[`j',15] = b_iv - b_lm 
		gen temp = ( e_iv*e_z - ( v_lm*z_lm + e_lm*e_z ) )/`xz_m'
		reg temp, robust
		matrix A[`j',16] = _se[_cons]
		drop temp	
		
		drop _all	
	}	
	//matrix list A
	
	clear
	svmat A, names(col)
	save ../result/sim`c'.dta, replace
	gen t1 = gap/gap_se
	gen t2 = lm/lm_se
	gen t3 = me/me_se
	forvalues j=1/3{
		gen r`j' = abs(t`j')>invnormal(0.975)
	}
	tabstat _all, col(stat) stat(mean sd p1 p5 p10 q p90 p95 p99) format(%8.4f)
}

